{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=workshops/ICRW_2023.ipynb)\n",
    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/workshops/ICRW_2023.ipynb)\n",
    "\n",
    "**An Introduction to Watershed Analysis with Leafmap and WhiteboxTools**\n",
    "\n",
    "This [notebook](https://leafmap.org/workshops/ICRW_2023) provides an introduction to watershed analysis with Leafmap and WhiteboxTools. It is designed for The Interagency Conference on Research in the Watersheds ([ICRW](https://icrwatersheds.org/)) 2023 workshop - [Working with Geospatial Hydrologic Data for Watershed Analyses in R and Python Using Web Services](https://icrwatersheds.org/icrw8/icrw8-schedule/june-5-workshops/).\n",
    "\n",
    "- Leafmap: <https://leafmap.org>\n",
    "- WhiteboxTools: <https://www.whiteboxgeo.com>\n",
    "- WhiteboxTools User Manual: <https://www.whiteboxgeo.com/manual/wbt_book>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "Uncomment and run the following cell to install necessary packages for this notebook, including leafmap, geopandas, localtileserver, rio-cogeo, pynhd, py3dep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install leafmap[raster] geopandas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import leafmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Create interactive maps\n",
    "\n",
    "Specify the map center, zoom level, and height."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[40, -100], zoom=4, height=\"600px\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Add basemaps\n",
    "\n",
    "Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_basemap(\"OpenTopoMap\")\n",
    "m.add_basemap(\"USGS 3DEP Elevation\")\n",
    "m.add_basemap(\"USGS Hydrography\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "Add NLCD land cover map and legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[40, -100], zoom=4)\n",
    "m.add_basemap(\"HYBRID\")\n",
    "m.add_basemap(\"NLCD 2019 CONUS Land Cover\")\n",
    "m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "Add WMS layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[40, -100], zoom=4)\n",
    "m.add_basemap(\"Esri.WorldImagery\")\n",
    "url = \"https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?\"\n",
    "m.add_wms_layer(\n",
    "    url,\n",
    "    layers=\"NLCD_2019_Land_Cover_L48\",\n",
    "    name=\"NLCD 2019 CONUS Land Cover\",\n",
    "    format=\"image/png\",\n",
    "    transparent=True,\n",
    ")\n",
    "m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## Get watershed data\n",
    "\n",
    "Let's download watershed data for the Calapooia River basin in Oregon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource=\"comid\", simplified=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Plot the watershed boundary on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_gdf(gdf, layer_name=\"Catchment\", info_mode=None)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "Save the watershed boundary to a GeoJSON or shapefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.to_file(\"basin.geojson\", driver=\"GeoJSON\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.to_file(\"basin.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## Download DEM\n",
    "\n",
    "Download a digital elevation model (DEM) for the watershed from the USGS 3DEP Elevation service. Convert the DEM to a Cloud Optimized GeoTIFF (COG)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.get_3dep_dem(\n",
    "    gdf, resolution=30, output=\"dem.tif\", dst_crs=\"EPSG:3857\", to_cog=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "Display the DEM on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"dem.tif\", palette=\"terrain\", layer_name=\"DEM\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "## Get DEM metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata = leafmap.image_metadata(\"dem.tif\")\n",
    "metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "Get a summary statistics of the DEM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata[\"bands\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "## Add colorbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_colormap(cmap=\"terrain\", vmin=\"60\", vmax=1500, label=\"Elevation (m)\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "## Initialize WhiteboxTools\n",
    "\n",
    "Initialize the WhiteboxTools class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt = leafmap.WhiteboxTools()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "Check the WhiteboxTools version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.version()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "Display the WhiteboxTools interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.whiteboxgui()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "## Set working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.set_working_dir(os.getcwd())\n",
    "wbt.verbose = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "## Smooth DEM\n",
    "\n",
    "All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.feature_preserving_smoothing(\"dem.tif\", \"smoothed.tif\", filter=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "Display the smoothed DEM and watershed boundary on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(\"smoothed.tif\", palette=\"terrain\", layer_name=\"Smoothed DEM\")\n",
    "m.add_geojson(\"basin.geojson\", layer_name=\"Watershed\", info_mode=None)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "## Create hillshade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.hillshade(\"smoothed.tif\", \"hillshade.tif\", azimuth=315, altitude=35)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "Overlay the hillshade on the smoothed DEM with transparency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"hillshade.tif\", layer_name=\"Hillshade\")\n",
    "m.layers[-1].opacity = 0.6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46",
   "metadata": {},
   "source": [
    "## Find no-flow cells\n",
    "\n",
    "Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.find_no_flow_cells(\"smoothed.tif\", \"noflow.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48",
   "metadata": {},
   "source": [
    "Display the no-flow cells on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"noflow.tif\", layer_name=\"No Flow Cells\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50",
   "metadata": {},
   "source": [
    "## Fill depressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.fill_depressions(\"smoothed.tif\", \"filled.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52",
   "metadata": {},
   "source": [
    "Alternatively, you can use depression breaching to fill the depressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.breach_depressions(\"smoothed.tif\", \"breached.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.find_no_flow_cells(\"breached.tif\", \"noflow2.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"noflow2.tif\", layer_name=\"No Flow Cells after Breaching\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56",
   "metadata": {},
   "source": [
    "## Delineate flow direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.d8_pointer(\"breached.tif\", \"flow_direction.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58",
   "metadata": {},
   "source": [
    "## Calculate flow accumulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.d8_flow_accumulation(\"breached.tif\", \"flow_accum.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"flow_accum.tif\", layer_name=\"Flow Accumulation\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61",
   "metadata": {},
   "source": [
    "## Extract streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.extract_streams(\"flow_accum.tif\", \"streams.tif\", threshold=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"streams.tif\", layer_name=\"Streams\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64",
   "metadata": {},
   "source": [
    "## Calculate distance to outlet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.distance_to_outlet(\n",
    "    \"flow_direction.tif\", streams=\"streams.tif\", output=\"distance_to_outlet.tif\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"distance_to_outlet.tif\", layer_name=\"Distance to Outlet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67",
   "metadata": {},
   "source": [
    "## Vectorize streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.raster_streams_to_vector(\n",
    "    \"streams.tif\", d8_pntr=\"flow_direction.tif\", output=\"streams.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69",
   "metadata": {},
   "source": [
    "The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use leafmap.vector_set_crs() to set the coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.vector_set_crs(source=\"streams.shp\", output=\"streams.shp\", crs=\"EPSG:3857\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"streams.shp\", layer_name=\"Streams Vector\", style={\"color\": \"#ff0000\", \"weight\": 3}\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72",
   "metadata": {},
   "source": [
    "## Delineate basins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.basins(\"flow_direction.tif\", \"basins.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"basins.tif\", layer_name=\"Basins\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75",
   "metadata": {},
   "source": [
    "## Delineate the longest flow path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.longest_flowpath(\n",
    "    dem=\"breached.tif\", basins=\"basins.tif\", output=\"longest_flowpath.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77",
   "metadata": {},
   "source": [
    "Select only the longest flow path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.select_largest(\n",
    "    \"longest_flowpath.shp\", column=\"LENGTH\", output=\"longest_flowpath.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"longest_flowpath.shp\",\n",
    "    layer_name=\"Longest Flowpath\",\n",
    "    style={\"color\": \"#ff0000\", \"weight\": 3},\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80",
   "metadata": {},
   "source": [
    "## Generate a pour point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81",
   "metadata": {},
   "outputs": [],
   "source": [
    "if m.user_roi is not None:\n",
    "    m.save_draw_features(\"pour_point.shp\", crs=\"EPSG:3857\")\n",
    "else:\n",
    "    coords = [-122.613559, 44.284383]\n",
    "    leafmap.coords_to_vector(coords, output=\"pour_point.shp\", crs=\"EPSG:3857\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82",
   "metadata": {},
   "source": [
    "## Snap pour point to stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.snap_pour_points(\n",
    "    \"pour_point.shp\", \"flow_accum.tif\", \"pour_point_snapped.shp\", snap_dist=300\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\"pour_point_snapped.shp\", layer_name=\"Pour Point\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85",
   "metadata": {},
   "source": [
    "## Delineate watershed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.watershed(\"flow_direction.tif\", \"pour_point_snapped.shp\", \"watershed.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"watershed.tif\", layer_name=\"Watershed\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88",
   "metadata": {},
   "source": [
    "## Convert watershed raster to vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.raster_to_vector_polygons(\"watershed.tif\", \"watershed.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"watershed.shp\",\n",
    "    layer_name=\"Watershed Vector\",\n",
    "    style={\"color\": \"#ffff00\", \"weight\": 3},\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
